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ABSTRACT 


A Two-Dimensional  Numerical  Model  of  Dry  Convection 
With  Three-Dimensional  Dynamics.  (August  1978) 
James  Charles  Weyman,  B.S.,  Grove  City  College 
M.S.B.A.,  Metropolitan  College 
Chairman  of  Advisory  Committee:  Dr.  Dusan  Djuric 


A numerical  model  for  the  study  of  dry,  three-dimensional, 
small  scale,  atmospheric  convection  is  presented  for  use  with  a 
two-dimensional  grid  in  the  vertical  xz  plane.  All  of  the  deriva- 
tives in  the  horizontal  y dimension  are  derived  through  the  assump- 
tions of  cyclostrophic  balance,  symmetry  of  circular  eddies,  and 
horizontal  isotropy  of  derivatives.  Variable  eddy  coefficients, 
proportional  to  the  deformation  field  and  the  square  of  the  grid 


interval,  are  used.  Comparisons  with  a similar  two-dimensional 
model  shows  the  thermals  in  this  pseudo  three-dimensional  model 
develop  sooner  and  reach  a larger  maximum  vertical  velocity  faster. 
These  results  are  comparable  to  what  other  researchers  of  three- 
dimensional  models  have  found.  The  economy  of  computation  of  this 
model  will  enable  other  investigators  to  make  better  use  of 
limited  computer  facilities. 


| 


limited  computer  facilities 
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1.  INTRODUCTION 

Numerical  modelling  is  one  of  the  promising  ways  to  study  the 
details  of  convective  flow  which  cannot  be  practically  observed.  The 
models  that  have  been  tested  thus  far  have  differed  greatly  in  the 
sophistication  of  the  cloud  physics  and  the  numerical  techniques 
employed.  One  of  the  earliest  experiments  was  carried  out  in  the 
mid-1950's  at  the  Los  Alamos  Scientific  Laboratory  at  the  instigation 
of  J.  von  Neumann.  The  results  of  this  model  were  later  published 
by  Blair  et  al.  (1959).  Blair's  model  involved  the  simulated  over- 
turning of  an  unstably  stratified,  two-layer,  incompressible  fluid 
system.  Then  Lilly  (1962)  proposed  a two-dimensional  slab-symmetric 
model  of  buoyant  convection.  This  model,  which  qualitatively  and 
quantitatively  resembled  the  convective  thermals  described  by 
Scorer  and  Richards  (1959) , did  not  exhibit  the  shape  preserving 
stage  assumed  in  theoretical  treatments  and  found  by  laboratory 
experiments.  Lilly  attributed  this  to  the  neglect  of  the  effects  of 
the  eddies  in  the  third  dimension.  Although  Lilly  met  with  limited 
success,  his  list  of  twelve  areas  to  be  studied  further  served 
as  a guide  for  other  researchers  in  the  investigation  of  convective 
processes  by  numerical  simulation. 

Using  many  of  Lilly's  suggestions,  a number  of  different  models 
were  investigated  in  the  1960's.  Ogura  (1962)  used  an  axially  symmetric 
model  to  simulate  a buoyant  mass  of  fluid  embedded  in  an  ambient 

The  format  and  style  of  this  thesis  follow  those  of  the  Journal 
of  Atmospheric  Sciences. 


fluid  of  uniform  density.  The  results  from  Ogura's  model  exhibited 
the  shape  preserving  stage  that  Lilly's  experiment  could  not.  How- 
ever, Ogura's  model  was  unable  to  handle  wind  shear.  Ogura  and 
Charney  (1960)  developed  a two-dimensional  slab-symmetric  model  to 
simulate  a squall  line.  The  serious  disadvantage  in  their  studies 
was  that  the  resulting  downward  motion  began  very  close  to  the  upward 
thermal,  thereby  cutting  off  the  supply  of  warm  air  necessary 
for  the  upward  motion.  The  thermal  was  greatly  weakened  after  this 
occurred.  Squires  and  Turner  (1962)  used  a one-dimensional  model 
to  study  cumulonimbus  updrafts.  They  allowed  for  the  entraining 
of  environmental  air  by  assuming  the  inflow  velocity  was  proportional 
to  the  upward  velocity  of  the  plume.  They  also  allowed  for  the 
incorporation  of  latent  heat  of  freezing  by  assuming  the  proportion 
of  ice  to  total  condensed  material  varied  linearly  with  temperature. 
Squires'  and  Turner's  model  yielded  results  which  seemed  reasonably 
consistent  with  observations,  especially  in  regard  to  cloud  shape. 
However,  since  the  model  was  one  dimensional,  it  could  not  adequately 
represeut  the  horizontal  variations  of  conditions  within  the  cloud. 
Orville  (1964,  1965)  used  a two-dimensional  slab-symmetric  model 
to  simulate  mountain  upslope  winds.  Although  his  results  were 
similar  to  the  observations  available,  he  believed  that  the  two- 
rather  than  the  three-dimensional  treatment  may  be  the  largest 
drawback.  The  reason  was  that  the  downward  vertical  motion  and  the 
energy  budget  were  not  sufficiently  satisfied  in  the  two-dimensional 
model.  Orville  (1968)  used  an  improved  two-dimensional  slab-symmetric 


model  to  simulate  the  development  of  cumulus  clouds  over  a mountain. 
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Although  Orville's  results  were  an  improvement  over  his  1964  and  1965 
work,  he  still  specified  his  major  problem  was  the  missing  third 
dimension. 

With  the  advent  of  larger  and  faster  computers  during  the  late 
1960's  and  early  1970' s,  it  became  possible  to  construct  three-dimen- 
sional models.  Deardorff  (1970)  was  one  of  the  first  to  model  in 
three  dimensions  when  he  investigated  three-dimensional  turbulent 
channel  flow  at  large  Reynolds  numbers.  Although  the  resolution 
was  modest  at  first  (24  X 14  X 20  grid  points) , the  cascade  of 
energy  from  large  scales  to  smaller  scales,  which  scarcely  occurred 
in  two  dimensions,  was  very  discernable  in  three  dimensions.  This 
turbulence  model  was  extended  by  Deardorff  (1972)  to  the  study  of 
neutral  and  unstable  planetary  boundary  layers  using  a grid  of 
40  X 40  X 20  points.  This  model  was  further  increased  to  40  X 40  X 40 
points  by  Deardorff  (1973)  to  study  the  use  of  subgrid  transport 
equations  in  a three-dimensional  model  of  atmospheric  turbulence. 
Following  Deardorff 's  work,  a number  of  other  investigators  began 
using  three-dimensional  models.  Fox  (1972)  used  a 12  X 12  X 54  grid 
in  a three-dimensional  model  to  simulate  a three-dimensional  thermal. 
Although  Fox  felt  that  three-dimensional  simulation  was  a great  improve- 
ment over  two-dimensional  work,  he  stated  that  turbulent  thermals 
could  not  be  completely  simulated  without  more  powerful  computers  and 
vastly  more  computer  time.  Other  three-dimensional  models  have 
been  proposed  by  Steiner  (1973),  Wilhelmson  (1974),  Miller  and  Pearce 
(1974) , and  Schlesinger  (1975) . Schlesinger  stated  that  the  number 
of  three-dimensional  models  has  not  been  larger  mainly  because  they 
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demand  large  amounts  of  computing  storage  and  expenditures  if 
adequate  resolution  is  to  be  achieved. 

From  the  above  discussion,  disadvantages  of  the  present  one-, 
two-,  and  three-dimensional  models  have  been  shown  to  limit  the 
complete  study  of  convective  processes.  An  alternative  model  was 
proposed  by  Deardorff  (1965) . He  described  a pseudo  three-dimensional 
model  which  allowed  the  third  dimension  to  be  simulated  while  the 
computations  were  done  in  two  dimensions.  This  method  would  combine 
the  advantages  of  the  two-  and  three-dimensional  models.  Due  to  poor 
results  when  treating  convection  at  small  Prandtl  numbers  and  to 
the  significant  advancement  in  computers  during  the  late  1960's 
(Deardorff,  personal  communication) , Deardorff  never  continued  this 
research,  but  instead  began  hxs  three-dimensional  model.  However, 
Gruneberg  (1975)  decided  to  continue  this  research  of  a pseudo 
three-dimensional  model  because  of  the  limitations,  suggested  by 
Schlesinger  (1975)  and  Fox  (1972),  of  three-dimensional  models  to 
most  researchers.  Although  Gruneberg  used  Deardorff' s idea  of  a 
pseudo  three-dimensional  model,  his  physical  and  statistical 
assumptions  for  the  simulated  third  dimension  (y)  were  different. 
Gruneberg 's  model  produced  results  which  were  only  slightly  different 
than  conventional  two-dimensional  models.  He  contributed  these 
results  to  his  treatment  of  the  pressure  force.  In  Gruneberg' s model, 
the  first  derivative  of  pressure  with  respect  to  y was  a statistical 
parameter,  and  not  based  on  the  properties  of  the  flow.  The  second 
derivative  of  pressure  with  respect  to  y was  not  included.  Jenkins 
(1976)  working  with  Gruneberg' s model,  remedied  this  situation. 
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The  statement  for  the  first  derivative  of  pressure  with  respect  to 
y was  based  on  the  assumption  that  for  convective  processes,  the 
magnitude  of  the  local  time  derivative  is  much  smaller  than  the 
centrifugal  force.  Therefore,  the  cyclostrophic  approximation  was 
used.  The  inclusion  of  the  second  derivative  of  pressure  was  based  on 
the  assumption  that  the  pressure  distribution  in  axially  symmetrical 
vortices  was  also  axially  symmetric. 

Although  Jenkins'  results  were  an  improvement  over  Gruneberg's 
findings  when  compared  to  three-dimensional  models,  serious  problems 
still  existed.  Numerical  instability  occurred  after  15  min  of  simu- 
lated time  which  severely  limited  the  usefulness  of  the  results. 

Also,  large  unexplained  gradients  in  the  vertical  motion  and  potential 
temperature  fields  developed  near  the  bottom  of  the  layer  under 
consideration.  At  present,  no  one  has  investigated  these  limiting 
features  of  this  potentially  valuable  model. 

It  is  the  purpose  of  this  research,  therefore,  to  determine 
if  the  two-dimensional  assumptions  may  be  replaced  by  pseudo 
three-dimensional  assumptions  in  a numerical  model  to  produce  thermal 
convection  in  better  agreement  with  true  three-dimensional  models. 
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2.  GOVERNING  EQUATIONS 

The  basic  equations  to  be  used  in  this  model  are  the 
horizontal  and  vertical  equations  of  motion,  the  thermal  diffusion 
equation,  and  the  continuity  equation.  These  equations  are  for  a dry, 
incompressible  atmosphere. 
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0 = u +v  +w 
x y z 


The  variables  u,  v and  w are  the  components  of  motion  in  the  x,  y 

and  z directions,  respectively.  The  use  of  subscripts  refers  to  the 

partial  derivative  of  the  subscripted  variable  with  respect  to 

the  subscript.  The  variable  p is  the  ratio  of  the  pressure  deviation 

from  a hydrostatic  basic  state  to  the  density,  b is  buoyancy,  F^,  F^ 

and  F,  are  the  friction  terms  and  F.  is  the  thermal  diffusion  term. 

3 4 

The  variables  u,  v,  w,  b and  p are  averaged  values  over  an  elemental 
grid  volume.  The  friction  terms,  the  thermal  diffusion  terms  and 
buoyancy  are  defined  as  follows: 


f\  = V • (K  Vu) 
x m 

F-  = v • (K  Vv) 
i m 

F,  = V . (K  Vw) 
J m 

F4  = V • (KjVb) 

b -g(-f) 
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The  variable  g is  the  acceleration  due  to  gravity,  0 is  a constant 
basic  potential  temperature,  and  is  the  deviation  from  the  horizon- 
tally (x  direction)  averaged  potential  temperature.  The  explanation 
of  this  choice  for  the  formulation  of  buoyancy  will  be  presented 
in  Section  3.  and  are  the  eddy  viscosity  and  thermal  diffusivity 
coefficients,  respectively.  These  coefficients  are  allowed  to 


vary  with  time  and  space,  and  are  computed  locally  by 

_ "bu.  7>u.  Du.  "bu. 

K = (cA)  1/2  (— -±-  + L-)  (— — i-  — -1-) 

m T>x.  ?>x.  T>x.  Tax. 

Di 


= 3. OK 


The  quantity  A is  the  representative  grid  interval,  c is  a dimension- 
less constant,  and  the  indexes  i and  j have  values  of  1,  2 and  3. 


This  formulation,  although  at  first  was  used  for  general  circulation 
models  (Leith,  1965,  Mintz,  1965,  and  others  ) , was  shown  to  be 
applicable  when  the  inertial  subrange  exists  on  scales  encompassing 
the  grid  interval  by  Lilly  (1967) , and  Leith  (1968) , and  Deardorff 
(1970).  The  formulation  of  and  will  be  discussed  more  fully 
in  Section  3. 

In  order  to  evaluate  the  pressure  field,  a balance  equation  for 
pressure  is  formed.  This  is  done  by  taking  the  divergence  of  ( 1) — (3) 


and  by  making  use  of  (5) . The  result  is 


p + p + p = div  r-A(V)  +F  +bk*l 

xx  yy  zz  ~ ~ 


' 


where  A (V)  represents  the  advective  terms  in  ( 1) — (3)  and  £ is  the  ver- 


tical  unit  vector.  Eq.  (13)  is  a Poisson  equation  of  the  form 


which  can  be  solved  directly  by  the  use  of  a Fourier  transform,  in 
which  the  two-dimension  problem  can  be  reduced  to  a set  of  one-dimen- 


sional problems. 
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3.  PSEUDO  THREE-DIMENSIONAL  ASSUMPTIONS 


In  this  model  all  of  the  computations  involving  the  various  terms 
are  done  in  the  vertical  x-z  plane.  Therefore,  all  of  the  y deriva- 
tives in  the  basic  equations  must  be  determined  using  other  equations 
or  model  assumptions.  Eq.  (5) , the  continuity  equation  for  an 
incompressible  atmosphere  was  used  once  in  the  derivation  of  (13) , 
the  balance  equation  for  pressure,  but  it  may  be  used  once  more 
in  order  to  eliminate  v . When  (5)  is  combined  with  (2) , the 

y 

resulting  equation  is 


v = -uv  +v(u  +w  ) -wv  -p  +F_  (14) 

t x x£  z y 2 

which  is  used  in  the  basic  equations  instead  of  (2) . 

The  second  derivative  of  v with  respect  to  y is  also  needed  in 
the  basic  equations.  When  the  divergence  of  the  advection  terms 
in  (13)  is  taken,  one  obtains  a term 


- (w  ) =-v  v - w . 

y y y y yy 


The  continuity  equation  is  used  again  to  eliminate  v . If  the 
partial  derivative  of  (5)  with  respect  to  y is  taken,  the  result  is 


v = - (u  ) - (w  ) . 

yy  X y z y 


Then  if  the  various  functions  are  continuous,  the  order  of  differen- 
tiation may  be  reversed  and  one  arrives  at 


v = - (u  ) - (w  ) 

yy  y x y z 


This  can  now  be  used  in  (13) 
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The  derivatives  u , w , and  b are  treated  as  new  variables  for 
Y y Y 

which  new  equations  are  needed.  These  new  equations  are  obtained 
by  taking  the  derivatives  of  (1) , (3) , and  (4)  with  respect  to 
y.  Again  assuming  that  the  various  functions  are  continuous,  the  or- 
der of  differentiation  is  reversed.  The  resulting  equations  cure 
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In  order  to  find  the  underlined  y derivatives  in  these  equations,  the 
assumptions  of  horizontal  homogeneity  and  horizontal  isotropy  in  free 
convection  are  made.  These  assumptions  mean  that  the  turbulence  has 
quantitatively  the  same  structure  in  the  two  horizontal  directions  of 
the  flow  field  (horizontal  homogeneity)  and  that  the  average  value  of 
any  function  of  the  velocity  components,  defined  in  relation 
to  a given  set  of  horizontal  axes,  is  unaltered  if  the  horizontal 
axes  of  reference  are  rotated  in  any  manner  (horizontal  isotropy) . 
Therefore,  by  making  the  assumption  that  the  flow  is  isotropic 
one  obtains 


-x  x 


-x  X 


2 2 
(U  ) = (V  ) 

y x 


2 2 
(w  ) = (w  ) 

x y 


and 


-x  x 


2 2 

tb  ) = (b  r 

y x 


(16) 


(See  Taylor,  1935) 
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The  Glossary  of  Meteorology  (1959)  says  that  although  atmospheric 
turbulence  is  generally  non-isotropic,  isotropic  turbulence  forms  the 
basis  of  most  theoretical  analysis  of  turbulent  flow.  Hinze  (1975) 
states  that  a knowledge  of  the  characteristics  of  isotropic  turbu- 
lence, notwithstanding  its  hypothetical  character,  may  still  form  a 
fundamental  basis  for  tne  study  of  actual,  non-isotropic  turbulent 
flows.  He  continues  by  saying  that  many  features  of  isotropic  turbu- 
lence apply  to  phenomena  in  actual  turbulence.  In  numerical  work, 
Deardor ff  (1965)  uses  the  horizontally  homogeneous  and  isotropic 
turbulence  assumption  in  a similar  manner  to  what  has  been  done  in 
this  research.  Therefore,  these  two  assumptions  are  not  expected  to 
be  too  restrictive. 

After  the  values  of  u , w , and  b are  determined  initially 

y y y 

at  the  first  time  step  or  are  found  at  the  beginning  of  each 
succeeding  time  step  by  the  left-hand  side  of  (15a),  (15b),  and 
(15c),  these  values  are  adjusted  so  that  the  isotropic  conditions 
are  satisfied.  This  is  done  by  correcting  u , w , and  b in  the 

y y y 

following  manner: 


Uy (new) 


u 

y 


w . , = w 

y(new)  y 


1/2 


I(wx>2 


Li<v2J 

i<y2- 

ny2 


1/2 


1/2 


(17) 


(18) 


b , . = b 

y(new)  y 


(19) 
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The  summations  that  are  used  in  (17)  - (19)  are  over  a horizontal 
row  of  grid  points  in  the  x direction.  These  equations  alter  u^,  w^, 
and  b in  such  an  amount  that  (16)  is  satisfied. 

y 

However,  if  horizontal  isotropy  is  assumed,  Deardorff  (1965) 
suggests  that  there  should  be  constraints  upon  the  lateral  (y) 
advection  so  as  not  to  affect  the  mean  values.  This  is  due  to  the 
fact  that  the  failure  to  constrain  the  lateral  advection  terms  could 
cause  spurious  instability.  This  will  be  discussed  more  fully  later. 

The  reason  for  these  constraints  arises  in  the  following  manner. 
The  equation  of  continuity  which  has  been  multiplied  by  -b  and  the 
advective  terms  of  the  thermal  diffusion  equation  are 


respectively. 

form. 


0 = -b  (u  +v  +w  ) 
x y z 


b,  = -ub  -vb  -wb 
t x y z 


If  these  two  equations  are  added  together. 


(21) 

(22) 

the  flux 


b = -(ub)  -(vb)  -(wb)  , (23) 

t x y z 

is  obtained.  If  this  were  a true  three-dimensional  model,  then  the 
average  at  each  level  over  the  horizontal  plane  (x-y)  would  be 


(24) 
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l 

I 

I 


With  cyclic  boundary  conditions  in  the  x and  y directions, 


= o 


and 


Wbf*y  - 0 


Therefore, 


■b^y  - q^T*y 


In  this  research  all  of  the  variables  are  defined  only  on  the  x-z 
plane,  and  an  average  over  the  x-y  plane  is  not  possible.  However, 
an  average  of  (23)  can  be  taken  in  the  x direction.  This  equation  is 


b = -(ub)  -(vb)  -(wb) 

t x y z 


(25) 


Here 


(ub)  x = 0 , 


(26) 


because  cyclic  conditions  are  assumed  in  the  x direction.  Since  the 
isotropic  assumption  assumes  similar  statistical  properties  on  the 
average  in  the  x and  y directions,  it  is  necessary  therefore  to  make 


(vb)y  ~ = 0 


(27) 
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This  requirement  is  analogous  to  the  three-dimensional  case  and 
assures  that  nothing  different  is  going  on  in  the  y direction,  on 
the  average,  than  in  the  x direction.  Although  it  is  true  that 
within  a single  x-z  vertical  slice  of  real  three-dimensional  motion  one 
does  not  expect  (27)  to  hold  precisely  unless  the  slice  extends  very 
far  in  the  x direction,  this  condition  should  not  be  violated  signi- 
ficantly or  systematically,  since  (26)  is  exactly  zero,  it  is 
only  consistent  that  (27)  be  made  to  hold  as  exactly  as  possible. 

Nonlinear  computational  instability  arises  because  of  the  non- 
linear interaction  between  different  wave  modes  in  the  numerical 
evaluation  of  the  advection  terms.  Lilly  (1965)  points  out  that  to 
insure  stability  in  the  absence  of  net  boundary  fluxes  the  finite 
difference  form  of  the  advection  terms  should  conserve  linear  and 
quadratic  quantities.  There  are  no  net  fluxes  in  the  x and  y 
direction  for  the  true  three-dimensional  case  (24)  shown  before, 
because  cyclic  conditions  are  assumed  in  x and  y.  In  this  case 
the  mean  values  are  correctly  conserved  in  the  horizontal,  since 
the  averaged  horizontal  flux  terms  in  (24)  are  zero.  For  the  pseudo 
three-dimensional  model  used  in  this  research,  there  are  no  net 
fluxes  in  the  x direction  because  of  cyclic  conditions,  and  the  x 
direction  flux  terms  average  to  zero  conserving  the  mean  values. 

The  y direction  flux  terms  should  also  conserve  the  mean  values  due  to 
isotropic  considerations  and  in  analogy  to  the  three-dimensional  case. 
There  is  no  way  to  ensure  this  requirement  unless  (27)  is  satisfied. 
This  condition  also  avoids  the  spurious  instability  mentioned  by 


15 


Deardorff  (1965).  A similar  argument  can  be  made  for  the  quadratic 
terms,  and  this  will  be  covered  later. 

Therefore  from  (27) , the  constraint  upon  the  lateral  advection 
so  as  not  to  affect  the  mean  value  of  b is 

vb  * = -bv  X . (28a) 

y y 

Similar  constraints  upon  the  lateral  advection  so  as  not  to  affect  the 
mean  values  of  u,  w,  u , w , and  b can  be  developed. 

y y y 

This  is  accomplished  by  multiplying  the  continuity  equation  first 
by  u to  receive  the  first  equation,  then  by  w to  get  the  second,  by 
u for  the  third,  by  w for  the  fourth,  and  by  b for  the  fifth. 

y y y 

Then  the  first  equation  can  be  added  to  the  advective  terms  of  (1) , 
the  second  equation  can  be  added  to  the  advective  terms  of  (3), 
and  the  third,  fourth,  and  fifth  can  be  added  respectively  to  the 
advective  terms  of  (15a),  (15b)  and  (15c).  An  average  in  the  x 
direction  of  the  five  resulting  flux  form  equations  is  taken. 

Because  cyclic  conditions  are  assumed  in  the  x direction,  and  the 
isotropic  assumption  assumes  similar  statistical  properties  on  the 
average  in  the  x and  y directions,  one  obtains 


x x 

VU  = -U  V 

yy  y y 


vw 


yy 


= -W  V 

y y 


vb 


yy 


-b  v 

y y 


■x 


■x 


(28b) 
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Another  important  constraint  is  that  for  p^.  If  the  continuity 
equation  is  multiplied  by  -p  and  then  added  to  the  pressure  force 
terms  of  the  kinetic  energy  equation,  the  result  is  a flux  form 
equation.  A spatial  average  in  the  x direction  of  the  resulting 
equation  is  taken.  Due  to  isotropy,  the  constraint  for  p^  is  given 


vp  = -pv 

y y 


This  constraint  is  very  important  if  p^  in  (14)  is  to  properly 
simulate  the  transfer  of  kinetic  energy  (generated  by  buoyancy) 
from  the  vertical  component  into  the  horizontal  component  of  motion. 

In  addition  to  the  constraints  upon  the  lateral  advection  so 
as  not  to  affect  the  means,  Deardorff  (1965)  suggests  that  there 
should  also  be  constraints  upon  the  lateral  advection  so  as  not  to 
affect  the  mean  squared  or  quadratic  values.  Even  if  the  assumptions 
connected  with  the  lateral  advection  terms  leave  the  x mean  unchanged, 
these  advection  terms  could  spuriously  alter  the  mean  squared 
values.  The  formulation  of  these  constraints  is  done  in  the 
following  manner.  If  (21)  is  multiplied  by  b and  (22)  by  2b  and 
the  resulting  equations  are  added  together  and  averaged,  one  obtains 


2 2 2 
• (ub  ) - (vb  ) - (wb  ) 

x y z 


Since  the  first  term  on  the  right-hand  side  of  (29)  equals  zero 


exactly  due  to  cyclic  conditions  in  x,  the  second  term  should  also 
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L 


i 


equal  zero  exactly.  This  is  necessary  because  of  the  horizontally 
isotropic  and  homogeneity  assumptions  as  pointed  out  before  and  to 
avoid  the  spurious  instability  as  discussed  previously.  Therefore 
the  constraint  upon  the  lateral  advection  so  as  not  to  affect  the 
mean  of  b2  is 


— 
v(b  )y 


(30a) 


Constraints  upon  the  lateral  advection  so  as  not  to  affect  the 
2 2 2 2 2 

mean  of  u . w , (u  ) , (w  ) , and  (b  ) can  be  developed.  These 

y y y 

constraints  are 


v(u2) 


2 

-u  v 


X 


y 


V (U2) 

y 


v(w2) 

y 


2 

-W  V 

y 


■x 


V (w2) 

y 


(b2) 


y y 


- -(v 


(30b) 


The  values  of  u , w , and  b , given  in  the  first  time  step  by 

y y y 

the  initial  conditions  or  after  that  found  for  a new  step  by  the  left- 
hand  side  of  (15a),  (15b),  and  (15c),  were  first  adjusted  for  the 
isotropic  conditions  by  (17)  - (19) . Now  these  values  must  be 
corrected  for  the  constraints  in  (28a) , (28b) , (30a) , and  (30b) . 

This  is  done  in  the  following  manner. 

Since  v is  determined  from  the  continuity  equation,  v can  be 

y y 

used  to  calculate 
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-bv 


and 


-b2v 


Then,  b^  can  be  corrected  to  satisfy  the  equations 


vb 


rx  and 


-x 


-x 


-bv 


v(b2)  or  2vbb  = -b2v 

y y y 


(31) 


The  necessary  corrections  for  b^  will  be  proportional  to  v and  2vb, 


b = b + k v 

y(new)  y 1 


b = b + k 2vb  , 

y(new)  y 2 


respectively,  where  k^  and  are  the  correction  factors.  These 
corrections  are  chosen,  because  it  cannot  be  expected  that  these 
averages. 


(32) 


vb 


and 


2vbb 


will  tend  to  the  correct  values  unless  one  of  the  members  (b  ) is 

y 
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2vbb  , 2vb(b  +k„2vb)X 

y (new)  = y 2 


x 

2vbb 

y 


2 2 
4v  b 


-x 


An  overrelaxed  iteration  procedure  is  used  to  correct  b , so  that 

y 

a b . , is  found  that  will  satisfy  both  equations  in  (31) . 

y(new) 

Similar  procedures  can  be  programmed  to  satisfy  all  of  the  con- 
straints in  (28a),  (28b),  (30a),  and  (30b).  Now  that  u , w , and 

y y 

^ have  been  adjusted  for  the  isotropic  conditions  (17-19)  and 
corrected  for  the  constraints  in  (28a) , (28b) , (30a) , and  (30b) , 
they  may  be  used  in  the  right-hand  side  of  (15a) , (15b) , and  (15c) 
to  predict  new  values  in  time  of  u , w^,  and  b respectively. 

Also  in  (15a) , (15b) , and  (15c) , values  of  u „ w , and  b are 

yy  yy  yy 

needed.  As  a result  of  enforcing  the  constraints  of  (28a) , 

(28b) , (30a) , and  (30b) , one  is  provided  with  a method  to  determine 
these. 

Egs.  (11)  and  (12) , the  formulation  for  the  eddy  coefficients 
used  in  this  model,  were  first  proposed  by  Smagorinsky  (1963) . 

Tn  this  formulation  the  eddy  coefficients  are  assumed  proportional 
to  the  magnitude  of  the  velocity  deformation  field  and  to  the  square 
of  the  grid  interval.  Lilly  (1967)  showed  that  this  treatment  is 
consistent  with  the  existence  of  a three-dimensional  inertial  sub- 
range on  scales  comparable  to  and  less  than  the  grid  interval.  In 

the  same  paper,  Lilly  (1967)  also  estimated  the  value  of  c in  (11) 
-3/4 

to  be  0.23of  where  Of  is  the  approximate  value  of  the  Kolmogorov 
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inertial-subrange  constant.  Deardorff  (1971)  shows  for  a r=  1.41, 

an  averaged  value  obtained  from  Pond  ct  al.  (1963),  Lilly's 

formulation  give  c = 0.176.  However,  Lilly  (1966)  found  that  when 

the  deformation  field  is  obtained  from  finite  differences  across 

single  grid  intervals,  c increases  25%  giving  c = 0.22.  Deardorff 

(1970)  found  a similar  value,  0.21,  for  the  case  of  an  unstably 

stratified  planetary  boundary  layer.  Due  to  the  staggered  grid 

arrangement  in  this  research,  it  is  necessary  to  take  finite 

differences  across  two  grid  intervals  for  about  50%  of  the  terms 

in  the  deformation  equation.  Therefore  it  is  found  that  a value  for 

c of  0.25  or  13%  larger  than  Lilly's  value  achieved  the  best  results. 

Deardorff  (1972)  improved  this  formulation  when  he  discovered  that 

the  ratio  of  K,  /K  had  to  be  between  two  and  three  to  avoid  excessive 
h m 

intensities  in  the  temperature  spectrum  at  large  wave  numbers.  A 
value  of  three  for  this  ratio  is  used  in  this  research. 

In  this  treatment  of  the  eddy  coefficients  and  the  diffusion 
terms,  an  assumption  is  needed  for  the  y derivatives  so  that  they 
may  be  included  in  this  pseudo  three-dimensional  model.  For  example, 
in  (1)  F^  is  defined  as 


F 


1 


(K  u ) 


m x x 


+ (K  u ) +(K  u ) 

m y y m z z 


From  previous  assumptions  u^  can  be  determined,  but  the  quantity 

(K  u ) cannot.  For  this  term  it  is  assumed  that 
m y y 


(K  u ) 
m y y 


(K  u ) 
m y x 
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at  each  point.  Since  horizontal  homogeneity  is  assumed,  (K  u ) 

m y y 

should  vary  in  space  and  time  in  a manner  consistent  with  (K  u ) . 

m y x 

This  assumption  fulfills  this  requirement,  whereas  the  other  two 
possibilities 

(K  u ) =0  , (K  u ) = constant  in  y direction 

m y y my 


(K  u ) = K u , K = constant  in  y direction 

m y y m yy  m 


do  not. 


Another  y derivative  that  must  be  determined  is  the  first 
derivative  of  pressure  with  respect  to  y.  This  is  needed  in  (14) 
to  calculate  v^.  For  the  convection  process  under  consideration,  the 
magnitudes  of  the  local  time  derivative  and  the  coriolis  force  are 
much  smaller  than  the  centrifugal  force.  Hence,  the  cyclostrophic 
approximation,  where  the  centrifugal  force  balances  the  pressure 
force,  can  be  used  to  calculate  p^.  The  cyclostrophic  balance  can 
be  stated  in  the  form 


-*p  - S' 


where  C is  the  three-dimensional  streamline  curvature  vector,  and  V 

r*J 

is  the  magnitude  of  the  three-dimensional  wind  vector.  The  curvature 
vector  points  toward  the  center  of  curvature  with  a magnitude  of 
1/r.  It  is  helpful  at  this  point  to  introduce  a natural  coordinate 
system.  The  direction  of  the  coordinates,  s,  n,  and  z,  in  this 


system  are  defined  by  unit  vectors  t,  n,  and  k,  respectively.  The  unit 
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vector  t is  parallel  to  the  flow  at  each  point,  £ is  normal  to  the 
flow  and  is  directed  toward  the  center  of  curvature,  and  k is 
normal  to  the  flow  and  equals  t X n . In  this  system  the  velocity 

A t r* 

may  be  written  as 

V = Vt  or  t = V/V  . 

A>  /W 

The  rate  of  change  of  t following  the  motion  may  be  derived  from 
geometrical  consideration  in  a manner  similar  to  that  used  by  Holton 
(1972)  to  find  dV/dt. 

Recalling  that  j^t|  =1,  then 

= &s/R  = |6t,J/|t|  =|6t|  = |£v  A | and  ( |6v  /v|)/  is  = 1/R  (34) 

can  be  shown  to  be  true  from  Fig.  1. 


&V/V 


Fig.  1.  A small  segment  of  flow  in  a natural 

coordinate  system  (After  Holton,  1972) . 

E 1 
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By  noting  that  ot  is  directed  parallel  to  n in  the  limit  as 

tsJ  /V 

0,  the  equation 


*bs 


(V  /V) 


/R  = C 


(35) 


is  obtained  from  (34).  After  performing  the  necessary  differentia- 
tion, the  curvature  vector  is  given  by 


D 

1 

"bV 
— * 

"bV 

“bs 

(V  /V) 
r-' 

= 7 

(36) 


The  derivatives  along  the  streamline  direction,  s,  are  obtained  by 
the  use  of  the  chain  rule  and  by  the  use  of  the  following  relations. 


*»V 

“av; 

"bs 

t* 

* ^ 

r 

= s (V/V) 

~bv 

“bv 

— br 

rJ 

*br 

= V/V 

*bs 

=”br 

' "bs 

~bV 


In  these  relations,  r is  the  cartesian  position  vector, 


is  a dyadic  tensor. 

r r 


1 


c = 

a/ 


r 


Eq. 

(V/V) 

n 


(36)  now  becomes 

T>vi  r 


V (V/V) 

A-'  A* 


and 


"br 


or 


£-7  [y71 -sv]^  • 

To  resolve  (37)  into  a m<-'  o workable  form,  the  approach  of 
Riesener  (personal  communic  an)  is  shown  below.  It  cam  be  shown 


that 
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— v v • V = — V (V  • V)  - — W • V 
v ~ ~ ~ V ~ ~ v — ’ ~ 


and  it  follows  that 


— V V • V = — 7 (V2)  = 7 V 

v ~ ~ Z 2v  X>  Z 


If  (38)  is  substituted  into  (37) , 


C = — V • 7 V - \ V V * 7 V • V 
^ v2  L"  ~ ~ v2  " " ~ ~ 

2 

is  obtained.  Since  V*V/V  = 1,  the  first  right-hand  term  of 
be  multiplied  by  this  to  produce 

V * V 


[,• 


7 V)  - 

**  **>  .2 


— V V • 7 V • V 

yZ  «-<  /*> 


i 


or 


C = 


[(V-7V)V-VV.7V 

A//V/A/  /W  ^ -V  /W  I 


If  A(V)  is  defined  as  V-VV  = A(V),  where  A(V)  represents  the 

A*  A#  A*  ^ A#  ^ ^ *■ 

advection  terms,  (40)  becomes 


C 


= — r A (V)V  - V A (V)  • 

v4  [j-  ~ ~ ^ ~ 


With  the  aid  of  the  vector  identity,  (a  x b)  x c = (ba  - ab) • 

<N A /%/  /*J  ^ I 

equation  for  C is  found  from  (41)  such  that 


S-7 


i V x A (V)  I 


X V 


(38) 

(39) 

(39)  can 


(40) 

(41) 

c,  an 

(42) 


With  the  aid  of  (33)  and  (42) , the  value  of  p^  may  be  specified  as 


-p  = V c_ 

y 2 


= ' ^u(uA2~vA^)  + w(wA2~vA3)J 


(43) 


and  this  value  can  be  used  in  (14) . 

The  last  y derivative  need  for  the  basic  equations  is  p . The 
second  derivative  of  pressure  with  respect  to  y is  essential  for  the 
three-dimensional  dynamics,  since  only  a balance  equation  for  pres- 
sure with  this  term  included  can  yield  a realistic  spatial  distribu- 
tion of  pressure.  An  estimate  of  p can  be  obtained  by  the 

yy 

assumption  that  the  pressure  distribution  in  axially  symmetrical 
vortices  is  also  axially  symmetric.  Then  the  available  elements 
of  the  flow,  V and  W,  give  an  opportunity  to  find  p . In  solid 

y y 

symmetrical  vortices  the  curvature  part  of  vorticity,  VC,  is  much 
greater  than  the  deformation,  def  V,  and  the  ratio. 


q = 


VC 


VC  + def  V 

gives  a measure  of  the  part  of  the  flow  that  can  be  represented 
by  a solid  vortex.  This  ratio  is  approximately  one  when  solid, 
completely  symmetrical  vortices  prevail.  Here  def  V is  the  familiar 
three-dimensional  deformation  of  the  flow. 


V-  [a 


2 2 2 2 2 2l  1/2 

def  V = 1 2(u  +v  +w  ) +(u  +v  ) +(u  +w  ) +(v_+w  ) 1 

xyz  y x zx  2 y 


1 


and  VC  is  the  magnitude  of  the  velocity  vector  multiplied  by  the 
magnitude  of  the  three-dimensional  streamline  curvature  vector  given 
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in  (42).  On  the  basis  of  symmetry  considerations,  p is  estimated  as 

YY 


P = (P  + P ) 
yy  xx  zz 


1/2 


VC 


VC  + def  V 


cos  a£ 


(44) 


where  of  is  the  angle  between  the  computing  plane  and  the  axis  of 

rotation  of  the  vortex.  This  estimate  is  realistic  as  can  be  seen 

from  the  examples  shown  in  Fig.  2.  For  the  case  shown  in  Fig.  2(a), 

p equals  p , and  p equals  zero  since  the  cos  90°  equals  zero, 
xx  z z yy 

for  case  2(b),  p equals  zero,  and  p equals  p . In  Fig.  2(c),  the 
z z xx  y y 

general  case  is  given  where  p is  defined  by  (44) . In  (44)  the  square 

root  of  the  quantity  in  parentheses  indicates  p will  be  about  equal 

to  the  larger  of  p or  p Therefore,  the  sign  of  p will  be  made 
xx  zz  yy 

the  same  as  that  of  either  p^  or  p depending  upon  which  quantity 
has  the  largest  absolute  value.  The  cos  may  be  calculated  using 

2 -2 
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1 2.^2  -2 


v+f+r 

where  , and  f are  the  components  of  the  rotation  vector 


(45) 


V X C = 1/V  V X A (V)  . 


With  the  use  of  (44)  and  (45),  the  balance  equation  for  pressure,  (13) 
may  now  be  written  as 
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p +p 
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In  summary,  the  necessary  assumptions  so  that  the  y derivatives 
may  be  included  in  this  pseudo  three-dimensional  model  are  (14) -(18), 
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(28),  (30),  (43),  (44),  and  the  assumption  for  (K^u^)  . With  these 
assumptions,  the  basic  equations  can  be  solved  step  by  step  in  time 
to  obtain  u,  v,  w,  and  b. 

After  the  y derivatives  have  been  introduced  in  the  described 
way,  it  may  be  useful  to  give  an  explanation  of  the  particular 
choice  of  the  formulation  of  buoyancy  in  this  model.  The  reason 
is  that  here  the  continuity  equation  does  not  ensure  that  the  vertical 
mass  flux  in  the  updrafts  is  matched  by  the  downward  flux  between 
the  thermals.  The  continuity  equation  is  here  identically  satisfied 
and  the  continuity  of  mass  must  be  controlled  by  the  formulation  of 
buoyancy.  If  the  buoyancy  is  defined  as  a deviation  from  a constant 
value,  difficulties  appear,  as  was  experienced  in  preliminary  com- 
putations. In  such  a case  it  depends  on  the  choice  of  initial 
conditions  how  w develops.  In  these  experiments  an  unstable 
stratification  is  used  initially.  Then  the  model  automatically 
assumes  that  the  lower  regions  are  warmer  than  the  adjacent  regions 
in  the  y direction.  Also,  the  upper  regions  are  assumed  colder  than 
their  laterally  adjacent  regions  in  the  y direction.  The  consequence 
is  that  all  lower  regions  acquire  a positive  vertical  acceleration 
and  the  upper  regions  a negative  one.  The  flow  develops  the  physically 
uninteresting  situation  where  whole  horizontal  belts  within  the  plane 
of  computation  arise,  while  other  belts  sink. 

Therefore,  in  order  to  simulate  ascending  and  descending  thermals 
side  by  side,  the  deviation  from  the  horizontal  average  is  used  as 
buoyancy  in  (3) . In  this  way  one  actually  uses  the  classical  parcel 
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method,  where  the  buoyancy  is  proportional  to  the  difference  in 
temperature  between  the  parcel  and  its  horizontally  adjacent  en- 
vironment, not  the  total  environment  (See  Haltiner  and  Martin,  1957, 
Gqs.  5- 5, 5-6). 

One  final  assumption  used  in  this  research  is  that  the  values  of 
u,  w,  b,  Uy,  w^,  and  b^  in  the  computational  xz  plane  are  near  the 
maximum  or  minimum  values  with  respect  to  the  y direction.  When 
compared  to  Deardorff  (1965) , this  is  a less  restrictive  assumption. 

He  assumed  that  the  values  in  the  computational  plane  were  the 
maximum  or  minimum  values.  It  was  found  through  a number  of  test 
runs  that  if  the  full  y advection  into  the  computational  plane  is  used 
the  values  in  the  v field  increase  significantly  more  with  time  than 
the  values  in  the  u field.  These  large  v values  then  dominate  all 
the  other  variables  in  the  various  equations.  However,  if  the  mag- 
nitudes of  the  y advection  terms  in  (1),  (3),  (4),  and  (15a)-(15c) 
are  decreased  to  one-half  their  original  value,  this  does  not  occur. 
This  decrease  is  quite  reasonable  if  the  magnitude  of  .the  variable  is 
already  near  a maximum  or  minimum,  because  one  would  not  expect  a large 
amount  of  advection.  Therefore  if  the  advection  is  large  it  should 
be  decreased,  and  a decrease  of  one-half  was  found  to  be  most 


beneficial. 
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4.  NUMERICAL  CONSIDERATIONS 

The  numerical  calculation  of  convective  flow  is  a rather 
sensitive  process  since  an  elliptic  equation  with  von  Neumann 
boundary  conditions  must  be  solved  at  each  time  step.  To  handle 
this  an  interlaced  arrangement  of  grid  points  is  used.  The  field 
of  computation  is  divided  into  computational  boxes  with  the  grid 
points  of  pressure  in  the  middle  of  each  box  and  the  velocity  com- 
ponents in  the  middle  of  the  sides  normal  to  the  respective  velocity 
components.  Since  all  calculations  are  done  in  the  xz  plane,  this 
scheme  leaves  v in  the  same  points  as  pressure.  The  buoyancy  is 
given  in  the  same  points  as  w so  that  vertical  accelerations  can  be 
evaluated  with  the  simplest  centered  differences.  The  sides  of  the 
computational  boxes  are  Ax  = Az  = 40  m,  so  that  the  distance  between 
the  nearest  grid  points  is  20  m and  the  distance  between  the  nearest 
grid  points  with  the  same  variable  is  40  m.  Fig.  3 shows  this 
staggered  grid  arrangement.  All  results  presented  here  are  done 
in  a field  of  20  x 20  computational  boxes. 

The  time  extrapolation  is  done  using  the  Adams-Bashforth  scheme, 
which  is 

The  time  step  is  variable  and  is  chosen  automatically  each  step.  The 
length  of  the  time  step  varies  accordingly  to  the  Courant-Friedrichs- 
Lewy  criterion  for  the  advection  terms  or  the  von  Neumann  criterion 
derived  for  the  diffusion  terms.  These  criteria  are  given  by 
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respectively,  where  A=  t±y,  = Az- 

Spatial  derivatives  are  approximated  by  upstream  differences 
for  advective  terms  and  centered  differences  for  diffusion  terms. 

When  centered  differences  are  used  for  the  advective  terms,  as  Jenkins 
(1976)  did,  large  alternating  temperature  gradients  are  found  near 
the  bottom  of  the  model  which  were  caused  by  the  colder  air  in  the 
descending  currents  and  the  constant  temperature  surface.  These 
gradients  persist  due  to  the  centered  differencing  scheme.  When  the 
upstream  differencing  scheme  is  used,  these  alternating  gradients 
do  not  appear. 

As  an  example,  the  finite  difference  formula  for  the  advection 
of  the  u-component  of  velocity  is  shown  here  as 

VU3 

( ) 

Ax  _ 


MV  " ' 


u + 2u  -u 
1 Q 3 


- [-V  <v] 
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vvvw6 


VU4  1 


(46) 


for  the  case  where  u and  w are  positive 


The  averaged  quantities 
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(v  +v,)/2  and  (w  +w„  +w,  +wd/4 

U ■>  Q 2 3 6 

are  specified  this  way  due  to  the  staggered  grid  arrangement,  whereas 
the  average  of 

(U1  +2UQ  +U3)/4 

is  necessary  to  prevent  non-linear  numerical  instability. 

The  lateral  boundary  conditions  are  periodic  in  x,  and  a solid, 
frictional  lower  boundary  and  a solid,  free-slip  top  boundary  are 
assumed.  The  boundary  conditions  for  the  lower  and  upper  boundaries 
are  given  by 

u = v = w = w = u = 0 at  z = 0 

y y 


and 


w=b  =u  = v = w =(u)  =0  at  z = z. 

z z z y y z top 

Harlow  and  Welch  (1965)  show  that  for  a frictional  boundary 
the  normal  velocity  at  the  first  grid  point  below  the  boundary  should 
be  equal  in  both  magnitude  and  direction  to  the  normal  velocity  at  the 
first  grid  point  above  the  boundary.  On  the  other  hand,  the  tangen- 
tial velocities  below  a frictional  boundary  should  be  equal  in 
magnitude  but  opposite  in  direction  to  the  tangential  components 
directly  above  the  boundary.  For  a free-slip  boundary  the  normal 
velocity  component  directly  above  the  boundary  should  be  equal  in 
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magnitude  but  opposite  in  direction  to  the  normal  component  directly 
below  the  boundary.  The  tangential  velocity  vectors  on  either  side 
of  a free-slip  boundary  should  be  equal.  These  requirements  are 
enforced  in  this  model. 

The  initial  fields  of  u,  v,  and  w are  set  equal  to  zero,  and  then 
they  are  overlaid  with  a perturbation  field  generated  by  a random 
number  program.  This  random  number  program  provides  a set  of  numbers 
for  each  level  whose  mean  is  zero  and  whose  range  is  ±0.1  ms1. 

The  initial  potential  temperature  field  is  obtained  by  setting 
the  surface  potential  temperature  to  280.2  K,  and  then  by  decreasing 
the  potential  temperature  at  a rate  of  2.5  K km  1 for  the  first 
seventeen  vertical  grid  intervals  or  to  a height  of  680  m.  Then 
from  680  m to  800  m,  the  potential  temperature  is  increased  at  a rate 
of  25  K ki«  1.  This  simulates  an  unstable  layer  beneath  a very  stable 
layer.  This  field  is  then  overlaid  with  a perturbation  field  with 
a mean  of  zero  and  a range  of  ±0.1  K for  each  layer.  Along  the  bottom 
layer,  which  represents  the  surface,  the  temperatures  are  not  per- 
mitted to  change.  This  procedure  allows  localized  hot  areas  at  the 
surface  to  persist  which  initiate  and  maintain  the  convective  thermals. 

For  u , w , T , u , w , and  T , where  T is  potential 

y y y yy  yy  yy 

temperature,  the  initial  fields  are  determined  in  the  following  manner. 
These  fields  are  first  set  equal  to  zero,  and  then  they  are  overlaid 
with  a perturbation  field  generated  by  a random  number  program.  This 


program  provides  a set  of  numbers  for  each  level  whose  mean  is  zero. 
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The  range  of  the  random  numbers  at  each  level  is 
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The  numerical  procedure  was  tested  for  numerical  stability 
using  a small  3X3  grid  network  over  an  extended  time  period. 

This  represents  a very  severe  test  since  there  are  eight  boundary 
points  and  only  one  interior  point.  The  test,  which  ran  for  1000 
time  steps  representing  slightly  over  two  hours  of  simulated  time, 
demonstrated  that  the  numerical  solution  is  stable. 
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5 . RESULTS 

The  computer  program  of  this  model  is  constructed  for  the  pseudo 
three-dimensional  model,  hereafter  referred  to  as  3D,  described 
earlier.  However,  if  one  variable  is  changed  in  this  program,  all 
of  the  y derivatives  are  set  equal  to  zero  and  the  result  is  a 
two-dimensional  slab- symmetric  model  denoted  here  by  2D.  This  is 
done  so  that  the  two  models  have  the  same  initial  and  boundary 
conditions.  This  allows  direct  comparisons  to  be  made  between  the 
2D  and  3D  cases.  These  comparisons  are  shown  in  Figs.  4-16. 

Figs.  4 and  5 show  the  vertical  velocity  and  potential 
temperature  deviation  fields  for  the  2D  case  after  approximately 
5 min  of  simulated  time.  Figs.  6 and  7 show  the  same  fields  for  the 
3D  case  at  a similar  time.  The  corresponding  diagrams  for  the 
2D  and  3D  models  are  very  much  alike.  However  upon  close  scrutiny, 
one  can  see  that  the  3D  case  (Fig.  6)  is  slightly  farther  along  in 
development  than  the  2D  case  in  Fig.  4 when  both  are  compared  to 
what  occurs  at  a later  time.  The  secondary  thermal  at  (19,6)  which 
will  dissipate  in  later  time  steps,  has  a smaller  vertical  velocity, 
0.238  * s \ in  the  3D  case  compared  to  0.418  m s 1 in  the  2D  model. 
The  downdraft  at  (11,5)  in  Fig.  6 which  will  also  dissipate  later, 
has  a minimum  vertical  velocity  of  -0.309  ms1,  compared  to  a value 
of  -0.359  m s 1 in  the  downdraft  at  (11,6)  in  Fig.  4.  The  maximum 
w in  the  whole  field  is  0.875  m s 1 in  the  3D  model  and  only  0.824  m 
s 1 for  the  2D  case.  At  about  10  min  of  simulated  time,  Figs.  8-11, 
the  difference  between  the  two  cases  is  more  apparent.  The  maximum 
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Fig.  4. 


2D  vertical  velocity  (10  m s x)  at  4 min  46  s.  Positive 
values  are  represented  by  solid  lines,  while  negative  values 
are  shown  by  dashed  lines.  The  tick  marks  and  numbers  in  this 
figure  and  the  succeeding  ones  along  the  left  hand  side  and 
the  bottom  show  the  position  of  the  vertical  and  horizontal 
grid  points,  respectively. 
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Fig.  11. 


3D  potential  temperature  deviation  ( 10  1 K)  from  280  K 
at  10  min  27  s. 
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w for  the  3D  model  is  17%  larger  and  80  m higher  than  that  for  the 
2D  model.  A secondary  thermal  at  horizontal  grid  point  4 in  Fig.  10 
is  better  organized  with  a maximum  w of  0.43  m s 1 in  the  3D  case 
compared  to  0.29  m s ^ in  the  2D  case.  The  downdrafts  are  also  far- 
ther along  in  development.  The  minimum  w in  the  downdraft  at 
horizontal  grid  point  18  in  Fig.  10  is  smaller  and  160  m lower 
than  that  for  the  2D  case.  The  major  difference  in  the  potential 
temperature  deviation  fields  at  this  time  between  the  two  models  is 
that  warmer  air  is  carried  to  a greater  height  in  the  3D  case.  In 
Fig.  11,  a value  of  -0.18  K reaches  the  600  m level  compared  to 
400  m level  for  the  2D  case.  Continuing  in  time,  Figs.  12-15  show 
the  two  fields  for  the  2D  and  3D  cases  at  about  15  min,  which  for  the 
2D  case  is  the  time  of  maximum  development.  For  the  2D  model 
there  are  one  updraft  and  two  downdrafts  giving  a transverse  roll 
circulation.  These  rolls  with  their  axes  normal  to  the  plane  of  the 
model  are  similar  to  the  circulation  patterns  that  Steiner  (1973) 
found  with  his  slab-symmetric  model.  He  stated  that  these  rolls 
would  not  necessarily  occur  if  three-dimensional  circulation  was 
possible.  This  is  confirmed  in  the  computed  3D  case,  since  the  roll 
circulation  is  greatly  modified.  The  vertical  motion  and  the  poten- 
tial temperature  deviation  in  Figs.  14  and  15  show  signs  of  being 
influenced  from  outside  the  computational  plane.  This  is  especially 
noticeable  for  the  two  updrafts,  one  at  (1,6)  and  another  at 
(6,12).  These  updrafts  were  initiated  by  warm  air  advection  into 
these  areas  from  the  simulated  y direction.  Likewise,  cold  advection 
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Fig.  15.  3D  potential  temperature  deviation  ( 10  1 K)  from  280  K 
at  15  min  24  s. 
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from  the  simulated  y direction  near  the  bottom  of  the  model  has 
enhanced  the  two  downdrafts  at  (5,1)  and  (18,1).  Another  note- 
worthy feature  of  Fig.  14  is  the  secondary  maximum  that  does  not 
appear  in  Fig.  12.  Soong  and  Ogura  (1973)  found  a similar  secondary 
maximum  in  their  axisymmetric  model  which  did  not  appear  in  the  slab- 
symmetric  model  used  for  comparison. 

The  time  variation  of  the  maximum  w values  at  each  time  step 
are  given  for  the  2D  and  3D  models  in  Fig.  16.  After  2 min  of 
simulated  time,  the  3D  maximum  w values  are  larger  than  those  of  the 
2D  model,  except  for  a brief  period  between  15  and  16  min.  The 

maximum  difference  occurs  at  12.7  min  when  the  2D  maximum  w is  only 

0.82  of  the  maximum  w in  the  3D  model.  In  order  to  compare  the 

two  models,  the  ratio  of  the  maximum  w for  the  2D  model  for  the  entire 

time  period  to  that  of  the  3D  model  is  formed.  In  this  research 
the  ratio  is  0.90,  while  similar  ratios  found  by  other  investigators 
who  have  compared  two-  and  three-dimensional  models  are  0.71 
by  Wilhelmson  (1972),  0.35  by  Steiner  (1973),  and  0.51  by  Wilhelmson 
(1974).  Although  these  values  are  quite  different,  they  affirm 
that  the  three-dimensional  maximum  w is  always  larger  than  the 
two-dimensional  case  in  similar  circumstances.  Fig.  16  also  shows 
that  the  maximum  w,  besides  being  larger,  occurs  sooner.  In  the 
3D  case  the  maximum  w occurs  at  12.7  min  with  secondary  maximums 
at  14.8  and  16.4  min.  In  contrast,  the  maximum  w in  the  2D  model 
is  not  reached  until  15  min  with  secondary  maximums  at  12  and  16  min. 


Therefore  the  2D  model  takes  18%  more  time  to  reach  its  maximum  them 


variation  of  maximum  vertical  velocity  for  3D  and  2D  models 
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the  3D  case.  This  is  in  good  agreement  with  other  researchers 
who  have  made  similar  comparison  of  this  time  difference.  They  have 
found  values  of  221  (Wilhelmson,  1974) , 7%  (Steiner,  1973) , and 
13%  (Soong  and  Ogura,  1973). 

Another  comparison  that  can  be  made  between  the  2D  and  3D 
models  is  the  rate  of  ascent  for  maximum  w.  The  level  of  maximum 
w for  the  2D  case  rises  at  0.50  m s ^ while  for  the  3D  case  the 
ascent  is  0.73  m s The  ratio  of  the  3D  case  to  that  of  the  2D 
is  1.46.  This  can  then  be  related  to  the  researchers'  work  mentioned 
before.  Soong  and  Ogura  (1973)  found  the  rate  of  ascent  for  the 
axisymmetric  to  be  2.24  m s \ while  the  rate  for  their  slab-symmetric 
model  was  1.52  m s The  ratio  of  these  two  is  1.47.  Wilhelmson 
(1974)  found  a similar  ratio  of  1.78.  The  lack  of  uniformity 
among  the  rates  of  ascent  among  the  various  researchers  could  be 
contributed  to  the  different  initial  or  boundary  conditions, 
whether  moisture  was  accounted  for  and  in  what  way  it  was  handled 
if  it  was  included,  or  different  lapse  rates.  However  the  important 
part  in  the  comparison  of  the  2D  and  3D  models  is  the  ratio  of  the 
ascent  rates,  and  the  ratio  found  in  this  research  compares 
favorably  with  others. 

From  these  comparisons,  one  can  see  that  the  3D  thermals  grow 
faster,  have  a larger  vertical  motion,  and  reach  this  larger  w value 
faster  than  the  thermals  in  the  2D  case.  Also  the  unlikely  trans- 
verse roll  circulation  formed  in  the  2D  case  is  greatly  modified 


in  the  3D  model 
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6.  CONCLUSION 

The  aim  of  this  study  was  to  show  the  possibility  of  including 
the  y derivatives  in  a two-dimensional  model  of  convection  in  such 
a way  that  a simulation  of  the  three-dimensional  flow  is  achieved. 
The  necessary  y derivatives  were  included  by  assuming: 


(a) 


(b) 

(c) 


horizontal  isotropy  of  derivatives  for  u , w , 
and  b , Y Y 

YY 

quasi-cyclostropic  balance  for  p , and 
symmetry  of  solid  vortices  for  p 


YY 


w , 
YY 


By  assuming  horizontal  isotropy  of  derivatives,  additional  constraints 
arise  so  that  the  simulated  y advection  does  not  affect  the  mean  or 
the  mean  squared  values.  Jenkins  (1976)  failed  to  constrain 
these  terms,  and  he  experienced  numerical  stability  problems.  After 
these  constraints  were  enforced  in  this  research,  the  numerical 
solution  was  stable. 

The  results  of  this  research  are  quite  encouraging  and  demonstra- 
ted that  the  two-dimensional  assumptions  may  be  improved  by  pseudo 
three-dimensional  assumptions  to  produce  thermal  convection  in  better 
agreement  with  true  three-dimensional  models.  However  certain 
limitations  need  further  investigation.  The  vertical  gradient 
of  the  vertical  velocity  is  stronger  above  the  velocity  maximum  than 
below  it.  This  is  probably  due  to  insufficient  mixing  in  the  model 
in  this  area.  Also  it  is  found  that  the  v value  directly  above  the 
maximum  w is  much  smaller  than  the  surrounding  values  which  restricts 
the  y advection  terms  there.  From  a number  of  other  experiments,  these 
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two  limitations  have  been  shown  to  cause  a slower  ascent  rate  and 
an  overall  smaller  maximum  w than  other  researchers  have  found. 

These  two  limitations  are  possibly  connected  since  v is  determined 

y 

by  the  continuity  equation  in  which  w^  is  a major  term. 

Another  limitation  is  that  no  lateral  spreading  of  the  thermals 
occur  as  they  rise.  Steiner  (1973)  reports  a similar  finding. 

He  felt  there  was  insufficient  dissipation  and  increased  the  value 
of  c which  is  used  in  the  calculation  of  the  eddy  viscosity 
coefficient  to  0.42.  This  value  or  one  between  the  present  value 
used  and  0.43  should  be  tested.  Experimental  research  is  also  needed 
to  determine  why  the  ratio  of  minimum  to  maximum  w in  this  model 
is  larger  than  those  values  reported  by  investigators  of  three-dimen- 
sional  models.  The  smaller  y advection  term  in  the  vicinity  of 
w maximum  mentioned  previously  may  be  the  cause  of  this. 

When  the  preceeding  limitations  are  resolved,  moisture  may 
then  be  added  and  the  model  extended  to  the  study  of  convection  which 
is  not  horizontally  isotropic.  This  extension  would  require  the 
introduction  of  empirical  constants  differing  slightly  from  unity 
which  will  relate  the  y derivatives  to  the  x derivatives. 

Should  the  simulation  of  atmospheric  convection  be  successful 
using  a pseudo  three-dimensional  model,  more  researchers  will  be 
able  to  investigate  this  process  more  fully  with  the  computer 
facilities  readily  available  to  them. 
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